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6.5 COMPUTATION OF DISCRETE QUADRATIC TFDs 0 

6.5.1 General Computational Procedure 

Article 6.1 deals with definitions and properties of the discrete WVD (DWVD) and 
other discrete quadratic TFDs. It shows that the general discrete quadratic TFD 
of an analytic signal z[n] is 

*M] = *£ E G\p, m] z[n—p+m] z*[n-p-m } e~^ km ' M (6.5.1) 

|m|<f 

= 2 DFT {G[n, to] * (z[n+m\ z*[n — m])} ; m£(M ) (6.5.2) 

where the support dimensions of the kernel do not exceed M samples in the lag (m) 
direction and P samples in the time (n) direction, and (M) means any set of M 
consecutive integers; cf. [1, p.444]. So the general procedure for evaluating such a 
TFD is: 

1. Formation of the instantaneous autocorrelation function (IAF) 

K z [n,m\ = z[n+m\ z*[n — m}; 

2. Discrete convolution in n (time) with the smoothing function G[n,m]\ 

3. Discrete Fourier transformation mapping m (lag) to k (frequency). 

For the DWVD, which has G[n,m] = 6[n}. step 2 reduces to an identity transfor- 
mation and may be omitted. The windowed DWVD has G[n,m\ = <5[n] g[m], so 
that step 2 reduces to multiplication of the IAF by g [to]. Some other quadratic 
TFDs, however, have special forms leading to computational procedures which are 
not degenerate cases of the above, and which may be simpler or faster. 

This article addresses some of the practical issues in computing quadratic TFDs 
of a real signal, examines various cases of the above procedure, and considers the 
spectrogram as one example of a special form leading to a simpler, faster algorithm. 

6.5.2 Computation of the Analytic Signal 

The usual definitions of quadratic TFDs, especially the WVD and the windowed 
WVD, assume an analytic signal in order to avoid interference terms between pos- 
itive and negative frequencies. For computational purposes, an analytic signal also 
avoids the need for 2x oversampling prior to computation of the IAF (see Sec- 
tion 6.1.1 and ref. [2]). So, given a real signal s(t), we must first compute the 
analytic signal z(t) associated with s{t). The simplest method is the direct ap- 
proach of filtering out the negative frequencies in the frequency domain. If a real 
signal s [n] is given for n = 0, 1, 2, . . . , TV — 1 and periodically extended with period 
TV, where TV is even (or is made even by zero-padding), the algorithm is: 
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1. Compute S[k\ = DFT{s[n]} for A: = 0,1,..., AT — 1 ; 

2. Compute ( S[jfe] for k = 0, f 

Z[k]=l 2 S[k] for jfc = l,2,...,f-l (6.5.3) 

[ 0 otherwise ; 

3. Compute z[n] = IDFT{Z[/c]} , where IDFT{. . .} denotes the inverse DFT. 

The treatment of the Nyquist term (k = N/2) and the precise meaning of “analytic” 
for a discrete-time periodic signal are explained in [3] ; these issues become significant 
if the signal has non-zero amplitude at the Nyquist frequency. Further details on 
implementation of TFDs, including computation of the analytic signal, are given 
in [4] . A time-domain algorithm for computing the analytic signal using FIR filters 
is described in [5]. 

6.5.3 Real-Time Computation of TFDs 

The formula for the discrete quadratic TFD [Eq. (6.5.2)] involves the expression 
z[n+m] where m is allowed to be positive, together with z*[n—m ] where m is 
allowed to be negative. The same applies to the DWVD 

W z [n, k] = 2 DFT {z[n+m] z *[n—m]} ; m£(N) (6.5.4) 

m—*k 

and the windowed DWVD 

W®[n, k\ = 2 DFT { g[m ] z[n+m] z*[n—m ]} ; me (M) (6.5.5) 

m — ► fc 

(both of these expressions are derived in Article 6.1). Both cases involve time- 
advanced signals; for any value of n, the computation of the TFD involves signal 
values up to z[n+ A], where A is some positive integer. In real-time computation, we 
cannot compute the TFD for time n until we know the signal values up to z[n-\- A]; 
thus A is the minimum latency of the computation. In the case of the DWVD 
[Eq. (6.5.4)], the range of m for which the IAF can be non-zero is maximized when 
n is at the center of the time-support of the signal; so the latency reaches a peak of 
half the signal duration. For the windowed DWVD, the latency is limited to half the 
window duration. For the general discrete quadratic TFD, the latency is limited to 
half the sum of the dimensions of the G matrix. The latency of the analytic signal 
computation must be added to that of the TFD computation. In all cases a smaller 
value of M not only reduces latency but also produces shorter FFTs, hence shorter 
computational delays; but the cost is reduced frequency resolution. 

Latency is one of two measures of merit for real-time computation of TFDs. 
The other measure is throughput, which depends on the efficiency of the numerical 
algorithms. Eq. (6.5.2) can be written 

p z [n, k\ = 2 DFT {R z [n, m]} 

m—*k 



(6.5.6) 
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where R z [n, m] = G[n, m\ * ( z[n+m\ z*[n—m\). Similarly, 

p z [n+l, k\ = 2 DFT {R z [n+1, m]}. (6.5.7) 

m — ► k 

The above two equations represent successive time-slices of the TFD. Multiplying 
the second equation by j, adding the result to the first equation and using the 
linearity of the DFT, we obtain 

p z [n, k] + jp z [n+ 1, k] = 2 DFT {R z [n, m] + jR z [n+ 1, m]}. (6.5.8) 

m — >k 



If the TFD is known a priori to be real, as it usually is, then Eq. (6.5.8) means that 
the successive time slices of the TFD are respectively the real and imaginary parts 
of the right-hand side, which involves only one FFT [6] . Thus the realness property 
can enhance efficiency by halving the required number of FFTs. It can also halve 
the storage requirement as it implies Hermitian symmetry in the smoothed IAF. 



6.5.4 Computational Approximations for Discrete-Time Kernels 

Table 6.5.1 reproduces the “G[n, to]” column of Table 6.1.2 (p. 240) and adds two 
special cases often found in the literature: BJ 1/2 denotes the Born-Jordan distribu- 
tion with a = 1/2 , while ZAM 2 denotes the Zhao- Atlas-Marks distribution with 
a = 2 . Some entries in the “G[n, m]” column of Table 6.5.1 call for continuous con- 
volution prior to sampling. At best, the evaluation of such a convolution in the 
time-lag domain requires oversampling. At worst, it requires the numerical evalu- 
ation of an improper integral arising from a singularity in G(t, r). In either case, 
computational inefficiencies will arise unless the smoothing effect of the convolution 
can be approximated in some other way. The chosen approximations, shown in the 
right-hand column of Table 6.5.1, are explained below. 

In the case of the B-distribution, the sole purpose of the convolution is to avoid 
aliasing. Without the convolution, and for typical values of the parameter (3 (e.g. 
/ 3 = 0.01), the time-lag kernel would be a continuous function with a narrow slot 
at m = 0 caused by the factor |2m|^. This factor is approximately unity for small 
nonzero values of to. The convolution fills in the slot, so that the factor is approx- 
imately unity at to = 0 also. This effect can be approximated by replacing |2 to| 
with [4m 2 + 1] 1 / 2 , as is done in Table 6.5.1. 

In the case of ZAM distribution, for a suitable (unbounded) to [to], the con- 
volution also ensures that G[n, 0] = 6[n ] , which in turn verifies the TM property. 
Without the convolution, we would have 



Gzam [n, to] = w[m] rect(f^) 



w[m\ if | an | < |2m| 
0 otherwise. 



(6.5.9) 



This gives G[n, 0] = tu[0] <S[n] , which verifies the TM property provided that 
to[0] = 1 . Accordingly, Eq. (6.5.9) is used in Table 6.5.1, although other approx- 
imations are possible. For example, we could sacrifice the TM property in favor of 
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Table 6.5.1: Computational approximations for time-lag kernels of selected discrete quadratic TFDs. 
In the “Distribution” column, subscripts indicate parameter values while the prefix “w-" means “win- 
dowed” by the function w[m\. For the spectrogram and tu-Levin distributions, the window w is as- 
sumed to be real and even. The “G[n,m\” column shows the exact kernels required for the avoidance 
of aliasing in the Doppler-frequency domain. If G[n,m\ cannot be computed as written, the “Approx.” 
column shows the suggested computational approximation. 



Distribution 

WVD 



G[n , m] 
<5[n] 



Approx. 



Levin 

BJ 



\8[n+m ] + i<5[n— m] 

[|4^ reCt (4^)] 

** [sine n sine m] 



/ |4a:ra|-|-l 

l o 



if |2n| < |4am| + 1 
otherwise. 



BJ1/2 



Modified B 
w-WVD 



[ T 2 mT rect ( 2 rrf ) ] 

** [ sine n sine m] 

cosh~ 2 ^ n 
^2 n cosh - n 

<5[n] w[m] 



w- Levin 



\w[m] (<5[n+m[ + 8[n—m]) 



ZAM 



[«’H rect (l^)] 

** [ sine n sine m] 



i 

1 2m | + 1 

0 



if |n| < |m| 
otherwise. 



w [m] if |an| < \2m\ 
0 otherwise. 



ZAM 2 



Rihaczek 

w-Rihaczek 

Page 

CW 

B 

spectrogram 



[«’H rect (2^)] 

** [ sine n sine m] 
8\n — to] 
w*[—m] S[n — m] 



8[n — \m\] 

\2m\ 4m 2 ) 

** [ sine n sine m] 



{ «)[»«] if |n| < m 
0 otherwise. 



5 M 

/ 7 T<J 

y 4m 2 +7r<T 




2 2 \ 
— 7T un \ 

4m 2 -(-7r(j J 



if m, = 0 
otherwise. 



|2m| l^ 3 

1 1 J, mno /wi 




yj 4m 2 + 1 


. 9 ^ OlilC III 

cosh 




cosh 2 n 



w[n+m] w[n—m ] 



some smoothing by using the approximation 

Gzam [n, m] « ^w[m] [l +tanh(|4m| — |2a?t|)] , (6.5.10) 



and we could salvage the TM property by using a separate definition for m = 0 . 
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For the Born-Jordan (BJ) and Choi-Williams (CW) distributions, the convolu- 
tions are needed to remove singularities at to = 0 and ensure that G[n, 0] = 5[n] . 
For the BJ distribution, we can remove the singularity and approximate the spread- 
ing in the [n, to] plane by replacing |4 qito| with |4am| + l . The result is 

G B j[n,m] « j 4 d^ iect (\S^+l) - ( 6 - 5 - n ) 

which is equivalent to the rule given in Table 6.5.1. For the CW distribution, a 
similar effect is obtained by replacing |2m| with [4m 2 + tot] 1 / 2 . This step, by itself, 
gives the kernel 

G cw [n,m] » exp(;=^^). (6.5.12) 

For n = to = 0 , this reduces to G[0, 0] = 1 , which is consistent with the re- 
quirement that G[n, 0] = 5[n] . However, for m = 0, Eq. (6.5.12) reduces to 
G[n, 0] = e -7m , which is only an approximation to <5[n]. Accordingly, a two-part 
definition of the kernel is used in the “Approx.” column of the table. With n = 0 , 
the kernel as defined in the table reduces to 

Gcw[0,to] = v /=4^I (6.5.13) 

which takes the value 1 at m = 0 and 1 / y/2 at m = ± y^cr/H . For realistic values 
of a (e.g. cr > 1), this gives a reasonable degree of smoothing in the m direction. 

An alternative approach to the problem of singularities, which is not pursued 
here, is to evaluate the kernels in the Doppler-lag [l,m] domain. This is efficient if 
we intend to evaluate the time-convolution by the FFT method, which also uses the 
[Z,to] domain. But it is still an approximation (if the time-lag form of the kernel 
is taken as the definition) because the analytical formulae for continuous FTs of 
standard signals are only approximations when applied to the DFT. 



6.5.5 Special Case: Direct Form of the Discrete Spectrogram 

The short-time Fourier transform (STFT) of the continuous-time signal x(t) with 
real window w(t) is defined (in Section 2.3.1) as 

F x(t,f) = F f {x(T)w(T~t)} = e~ j2nft JF f {x(T + t)w(T)} 

/ OO 

x(t + t) w(t) e~^ 2 ^ T dr. 

-OO 

It is shown in Chapter 2 that the spectrogram S™(t, /), which is simply the squared 
magnitude of the STFT, can also be considered as a quadratic TFD with kernel 
w(t+ 1 ) w(t— §) . The discrete form of this kernel is w[n+m] w[n — m] . Hence the 
discrete spectrogram can be conveniently evaluated using the general procedure 
described in Section 6.5.1 above. But it is simpler and more efficient to discretize 
the continuous spectrogram directly. 



(6.5.14) 

(6.5.15) 
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Theorem 6.5.1: If the spectrogram S™ is modified by ideally sampling w(t) at 

t = m/f s (6.5.16) 

where m is an integer and f s is the sampling rate, and if 

w(t) = 0 for |t| > (6.5.17) 

where M is a positive integer, and if the modified TFD is denoted by S™ , then 

I |2 



GW ( n kf B \ _ \ m ( m+n\ m\ „-j2irkm/M 

x y~> w ^> e 



(6.5.18) 



\m\<M/2 



Proof/explanation: When w(r) is sampled, the integrand in Eq. (6.5.15) becomes 

oo 

x(t + t) w(t) e ~ ■ ?27r ^ T 5 (t — j 1 ) (6.5.19) 

m= — oo 

so that the STFT becomes 

oo 

F™(t,f) = e - j2nft ]T + i) ) e" J ' 2,r/m//s . (6.5.20) 

m= — oo 

By Eqs. (6.5.16) and (6.5.17), the summation is restricted to \m\ < M/2 , giving a maximum of M 
terms. 1 The sampling in r makes F™(t, f) periodic in / with period f s , while the time-limiting 
in r gives a frequency resolution of M bins per period. So it is convenient to let 

/ = kf s /M (6.5.21) 

where k is an integer. With these restrictions, Eq. (6.5.20) becomes 

F”(t, ^) = e -:>^kf s t/M x(^ -ft) w(^) e -J^km/ M (6.5.22) 

\m\<M/2 

Putting t = n/f s to match the quantization of r, then taking the squared magnitude of the discrete 
STFT, we obtain Eq. (6.5.18). ■ 



With a change of notation, Eq. (6.5.18) becomes 



S™[n, k\ 



x[m+n\w[rn\e- j2wkm ' M 

m\<M/2 



2 



(6.5.23) 



This S™[n,k] is the discrete spectrogram of the discrete-time signal x[n] with 
window w[m\. If the summand is extended periodically in m with period M (i.e. 
extended periodically in r with period M// s ), we obtain 



S™[n,k\ 



x[rn+n\w[m\e~ j27rkm/M 

me(M) 



2 



(6.5.24) 



1 M terms for odd M; M — 1 terms for even M . 
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where (M) denotes any set of M consecutive integers . 2 This may be written 



S™[n,k]= DFT {x[m+n] w[m]} ; toG(M). 



(6.5.25) 



The time support of S™[n, k] is that of x[n] * w[n ] , corresponding to x(t) * w(t) . 
If this has a duration not exceeding N samples, then the non-zero elements of the 
discrete spectrogram may be represented by an N x M real matrix. Only half of 
the M columns are needed for the non-negative frequencies, which are sufficient if 
x(t) is real. 

Eq. (6.5.23) involves a;[n + m] where |m| < M/2 and M is the window length 
in samples. So, in real-time computations, the latency of the discrete spectrogram 
computed by this formula is half the window length. 



6.5.6 Sample Code Fragments 

In view of the current popularity of Matlab™, we illustrate this Article with 
some code fragments from the experimental Matlab function tlkern.m, which 
computed all of the TFDs plotted in Article 5.7. The input parameters of the 
function specify the kernel in terms of a time-dependent factor g\ [n ] , a lag-dependent 
factor g 2 [m\, and an “auxiliary factor” <73 [n, m}. The overall time-lag kernel G[n, m\ 
is then computed as 

G[n,m] = g 2 [m\(g 1 [n] *g 3 [n,m]) = (g 2 [m]g l [n}) *g 3 [n,m}. (6.5.26) 

This scheme allows the computation of a wide variety of quadratic TFDs in un- 
der 320 lines of code, including exceptions for direct computation of the discrete 
spectrogram. 



6. 5. 6.1 Example 1: MBD (General Algorithm) 

For a separable kernel, the auxiliary factor would normally be omitted (i.e. taken 
as S[n, to]), while the time-dependent and lag-dependent factors would have input 
parameters specifying their types (e.g. Hamming or Hanning) and durations (in 
samples). Although the kernel of the modified B-distribution (MBD) is separable 
(see Article 5.7 and Table 6.5.1), its parameter /? is not a duration. The MBD 
kernel is therefore specified using the factors 

9i[n] = <5M ; 32 M = 1 ; g 3 [n,m] = n • (6.5.27) 

Notice that the auxiliary factor is the complete kernel. 

To compute the MBD in Fig. 5.7.2(e) on p. 220, the function tlkern is called 
with the following significant parameters: 



2 For even M, the periodic extension is padded with a zero term. 
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s = signal vector 
N = 128 = assumed period 
tr = 1 = time resolution 
tf = 'delta' = string specifying g\ [n] 

If = 'l' = string specifying g 2 [m] 

af = 'mb' = string specifying form of g% [n, m] 

ap = 0.2 = auxiliary parameter (/?). 

All internal computations, including IAF generation, are designed to be valid for 
periodic signals. Therefore, to compute the IAF of a non-periodic signal such as the 
one in Fig. 5.7.2(e), the assumed period N must be at least twice the signal length 
to avoid wrap-around effects. Because the time support of the IAF is identical to 
that of the signal, the same value of N is also sufficient to avoid wrap-around in the 
subsequent convolution with g\ [n] . 

The output is the real matrix tf d(l : Mpad+1 , 1 :Nsel) , whose dimensions Mpad+1 
and Nsel are assigned early by the statements 

Mpad = 2~ceil(log(2*M)/log(2) ) ; 7« lag-to-f requency FFT length 

Ncut = min(N,length(s) ) ; % duration of TF plot 

Nsel = ceil(Ncut/tr) ; °/ 0 no. traces in TF plot 

where M is the support length of the kernel in the lag direction; in this case M has 
been set to length (s) because of the constant “lag-dependent” factor. 

Preliminaries: The analytic signal is computed by the frequency-domain method. 

If N is even, the Nyquist term has Matlab index N/2+1 and the amplitude at that 
frequency is left unchanged [3] . If N is odd, there is no Nyquist term. The following 
code handles both cases: 

Noff = fix(N/2); 

z = fft (real(s) ,N) ; */* s truncated or padded 

z(2:N-Noff) = 2*z(2:N-Noff ) ; */, positive frequencies 

z(Noff+2:N) =0; */, negative frequencies 

z = if ft (z) ; 

For this kernel, the time-dependent factor gl and the lag-dependent factor g2 are 
computed by the statements 

gl ( 1 :N) = 0; 

gl(l) = l; 

g2(l:Mpad) = 1; 

where “ . . . ” denotes one or more line(s) of control code, or code that is skipped in 
this case. The auxiliary factor g3 (the whole kernel in this case) is computed by 

Moff = fix(M/2); 

g3 ( 1 : N , 1 : Mpad) = 0; 
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temp(l:N) = 0; 
for n = -Noff:Noff 

temp(l+rem(N+n,N) ) = (cosh(n) ) ~ (-2*ap) ; 
end 

temp = temp/ sum (temp) ; °/ t normalize 
for m = -MoffrMoff 

g3( : , l+rem(Mpad+m,Mpad) ) = temp.’; 
end 

where ap denotes the auxiliary parameter (/?), and the remainder (rem) function 
causes high array indices to represent negative values of time and lag. 

Step 1 — Formation of the IAF: The IAF matrix K(1 :N, 1 :Mpad) is formed by 

for n = 1:N 

for m = -Moff:Moff 

7# K(n,m) = z(n+m)z'‘*(n-m) , with corrected indices: 

K(n, l+rem(Mpad+m,Mpad) ) = z(l+rem(2*N+n+m-l ,N) ) *conj (z(l+rem(2*N+n-m-l ,N) ) ) ; 
end 
end 

where the “corrected” indices allow handling of periodic signals. 

Step 2 — Convolution in time: The assembly of the time-lag kernel and the con- 

volution in time with the IAF are performed together. The smoothed IAF is 

R z [n,m] = K z [n,m\ *G[n,m ] = K z [n,m] * ( 32 [tn] gi[n] * < 73 ( 71 , m]). (6.5.28) 

The above convolutions may be taken as circular if the assumed period is sufficiently 
long, in which case 

{ DFT 

This is implemented by the following code, in which K( : ,mcorr) is initially the m th 
column of the IAF, but is overwritten by the m th column of the smoothed IAF: 

for m = -Moff:Moff 

mcorr = l+rem(Mpad+m,Mpad) ; 



R z [n,m\ = IDFT 



{K z [n,m]} DFT {g 3 [n,m\} DFT {gi[n] g 2 [m]}\ . (6.5.29) 



K( : , mcorr) = if f t (f f t (K( : , mcorr) ) . *f ft (g3( : , mcorr) ) . *f f t (gl . ’ *g2 (mcorr) ) ) ; 



end 

(The factor < 72 \rn\ could be taken outside the IDFT, but this would not improve the 
efficiency of the code because g2 (mcorr) is a scalar.) The FFT method of convo- 
lution is useful in experimental code because of its generality, but is not necessarily 
the most efficient method, especially if one of the convolved sequences is short. 
Now we apply the time-resolution (tr): 

for nsel = l:Nsel 

7, nselth column of r is selected row of K: 
n = l+tr*(nsel-l) ; 
r ( : ,nsel) = K(n, : ) . ’ ; 
end 
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Step 3— DFT: The final DFT (lag to frequency) is computed by 

r = f ft (r) ; 

which, for the sake of generality, does not take advantage of realness. 

Final adjustments: The following code scales the TFD and repeats its first row 

(the zero-frequency row) so that the TFD spans a full cycle in the frequency domain: 

tfd = [real (r) ;real(r(l ,:))]. *(Ncut/Nsel/Mpad) ; 

The scaling ensures that the sum of the matrix elements is close to the signal energy 
regardless of the time resolution. 

6. 5. 6. 2 Example 2: Spectrogram (Special Case) 

The spectrogram in Fig. 5.7.3(f) on p. 221 was computed by the same function tlk- 
ern. For the spectrogram, the parameters s, N and tr are the same as for the MBD, 
while tf is ignored. Other significant parameters are 

If = 'rect.' = string specifying type of window 
M = 17 = window length (in samples) 
af = ' = string calling for spectrogram. 

The output is tf d(l : Mpad/2+1 , 1 : Nsel) , where the dimensions are assigned as for 
the MBD, except that the window duration M is read as an input parameter and 
not overwritten. 

The analytic signal is computed as for the MBD, although this is not strictly 
necessary for the spectrogram. 

The rectangular window is computed by 

Moff = fix(M/2); 
g2(l:Mpad) = 0; 

for m = -Moff:Moff 

g2(l+rem(Mpad+m,Mpad) ) = 1; 
end 

The matrix K(1 :N, 1 :Mpad) normally represents the IAF, but for the spectrogram 
it is assigned differently: 

for n = 1:N 

for m = -Moff:Moff 

°/ 0 K(n,m) = z(n+m)g2(m) , with corrected indices: 

K(n, l+rem(Mpad+m,Mpad) ) = z(l+rem(2*N+n+m-l,N))*g2(l+rem(Mpad+m,Mpad)) ; 
end 
end 

The code that applies the time resolution and perforins the final DFT (lag to 
frequency) is the same as for the MBD. But the final adjustment is different: 

tfd = (abs(r(l : Mpad/2+1, :))). ''2.* (Ncut/Nsel/Mpad/sum (g2. ''2)) ; 
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The magnitude-squared operation alters the relationship between the window and 
the scaling of the TFD. Also note that the above step uses only half the columns of 
the Fourier-transformed r matrix, namely those corresponding to the non-negative 
frequencies. Efficiency could be further improved by exploiting the analytic signal 
to halve the sampling rate. 

6.5.7 The TFSA package 

The Time-Frequency Signal Analysis ( TFSA) package is a set of functions developed 
over more than a decade at the Signal Processing Research Centre, Queensland Uni- 
versity of Technology, for computing modulated signals, quadratic and polynomial 
TFDs, ambiguity functions, wavelet transforms and scalograms, and various esti- 
mates of instantaneous frequency. As this is a production package rather than an 
experimental package, computationally intensive functions are precompiled and op- 
timized for efficiency, and an interactive user interface is added. The current version 
is distributed as a Matlab toolbox, so that TFSA functions can be used with other 
computational and graphical functions of Matlab. Further information is available 
at http://www.sprc.ciut.edu.au/ or http://www.eese.bee.ciut.edu.au/research/spr/ . 

6.5.8 Summary and Conclusions 

High-level programming languages with built-in FFT functions and matrix opera- 
tions have made it possible to construct compact yet highly versatile functions for 
computing quadratic TFDs. Use of a common algorithm for all TFDs is convenient 
for the programmer. But, as illustrated by the direct form of the spectrogram, ef- 
ficiency can sometimes be improved by using different algorithms in special cases. 
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